Identification of unique genomic signatures in patients with fibromyalgia and chronic pain

Fibromyalgia (FM) is a chronic pain syndrome characterized by widespread pain. The pathophysiology of fibromyalgia is not clearly understood and there are no specific biomarkers available for accurate diagnosis. Here we define genomic signatures using high throughput RNA sequencing on 96 fibromyalgia and 93 control cases. Our findings revealed three major fibromyalgia-associated expression signatures. The first group included 43 patients with a signature enriched for gene expression associated with extracellular matrix and downregulation of RhoGDI signaling pathway. The second group included 30 patients and showed a profound reduction in the expression of inflammatory mediators with an increased expression of genes involved in the CLEAR signaling pathway. These results suggest defective tissue homeostasis associated with the extra-cellular matrix and cellular program that regulates lysosomal biogenesis and participates in macromolecule clearance in fibromyalgia. The third group of 17 FM patients showed overexpression of pathways that control acute inflammation and dysfunction of the global transcriptional process. The result of this study indicates that FM is a heterogeneous and complex disease. Further elucidation of these pathways will lead to the development of accurate diagnostic markers, and effective therapeutic options for fibromyalgia.

such as TGF-beta, IL-6, IL-21 and IL-23 that promote Th17 differentiation confirming the presence of chronic inflammatory process 17 .Proteomic analyses of plasma proteins from women with FM have shown increased levels of specific plasma proteins suggesting systemic differences in protein expression associated with different clinical parameters 18,19 .A hypomethylated DNA pattern enriched for genes implicated in stress response and DNA repair has also been reported in FM 20 .However, the findings of many of these studies are often contradictory, including a small number of samples and many have not been independently confirmed.Finally, Chronic pain may not always be linked to genes but may be caused through interaction between genes, environmental and life style factors 21 .Due to significant phenotypic heterogeneity, multi-omics approaches may shed light on the underlying biological processes involved in FM and chronic pain.
To gain a comprehensive overview of the transcriptional processes and to define a potential genomic signature of FM, we performed high throughput RNA sequencing (RNA-seq) on peripheral blood mononuclear cells (PBMC) from 96 FM patients and 93 control individuals.Inclusion in the FM group required meeting the criteria of the 2016 College of Rheumatology and having a positive FM cytokine assay 3,6 .Our major objectives were to identify transcriptional differences between FM and healthy controls; to identify subgroups within the FM cohort, and to improve diagnosis and patient stratification using FM specific biomarkers.First, we assessed differentially expressed genes in FM cases compared to healthy controls.We then used bioinformatic tools to identify subgroups of FM patients with distinct genetic signatures.

Study participants
This study was performed with the approval of the institutional review board of the University of Illinois at Chicago (Office for the Protection of Research Subjects, OPRS) and all methods were performed in accordance with the relevant guidelines and regulations.The study groups consist of 96 FM patients (91 females and 5 males) and 93 control cases (41 females and 52 males).The participants (FM and controls) included in this study did not overlap with the participants from the previous study published in 2012 6 .All participants provided written informed consent.The inclusion criteria for FM followed the 2016 criteria of the American College of Rheumatology [3][4][5] and had a positive fibromyalgia assay (FM/a) 3,6 .The control group did not fulfill the 2016 criteria of the College of Rheumatology and had a negative fibromyalgia assay (FM/a).The FM/a included expression analysis of four cytokines, IL6, IL8, Mip1-α/CCL3 and Mip1-β/CCL4 and relied on the functioning of viable PBMC as previously described 6 .Exclusion criteria, both for patients and controls, were the presence of any other chronic disease (diabetes, heart disease or cancer).A questionnaire was used for the collection of demographic and clinical data from participants.None of the patients were treated with anti-inflammatory drugs at the time and 3 months before the start of the study.

Sample collection, RNA extraction, library construction and sequencing
For genomic analyses, blood samples (9-10 mL) from FM and control individuals were collected in Streck tubes (Streck, La Vista, NE).Samples were centrifuged at 2500 rpm to separate the plasma, PBMC and RBC layers.After removing plasma, the RBCs were lysed using Qiagen RBC lysis buffer and centrifuged for 20 min at 2500 rpm at room temperature.Cell pellets were washed in PBS twice to remove any trace of RBC, then homogenized using a QIAshredder homogenizer in 600 μL RLT buffer containing β-mercaptoethanol.Total RNA was extracted using QIAamp RNA blood mini kit (Qiagen, Germantown, MD).RNA quality was assessed using Agilent RNA screen tape and TapeStation 4200 (Agilent Technologies, Santa Clara, CA) and a nanodrop spectrophotometer was used to estimate the concentration and purity of RNA (Nano-Drop Technologies, Wilmington, DE).
RNA libraries were prepared using Agilent SureSelectXT RNA Direct workflow following manufacturer's instructions.Briefly, 200 ng of RNA was transferred into strip tubes and samples were completely dried at 30 °C in a vacufuge (Eppendorf).RNA-seq fragmentation mix was added to each sample, mixed by vortexing gently at 2000 rpm for 10 s.Fragmentation was carried out at 94 °C for 8 min and kept at 4 °C.RNA-seq first strand master mix was added to each sample and the first strand was synthesized using the following thermal cycling condition 25 °C (10 min), 37 °C (40 min) then maintained at 4 °C with heated lid on.The first stand was purified using Agencourt AMPure XP beads (Beckman Coulter Genomics).Second-strand cDNA was synthesized and end-repaired at 16 °C for 1 h.The second-strand cDNA was purified using Agencourt AMPure beads.The 3' ends of cDNA were dA-tailed at 37 °C for 30 min and adaptors were ligated to each dA-tailed cDNA at 20 °C for 15 min.The adaptor ligated cDNA was purified using Agencourt AMPure beads and was amplified using pre-capture thermal cycling conditions.The quality of the pre-captured library was assessed using D1000 ScreenTape on Agilent 4200 TapeStation system.The region between 150 and 400 bp was used for quantification.Hybridization was carried out overnight at 65 °C in a thermal cycler using 200 ng of pre-captured library.For capturing the targets, SureSelectXT Human All Exon V6 baits were used (Agilent Technologies Inc., Santa Clara, CA).The captured libraries were then amplified to add the index tags and were purified using Agencourt AMPure beads and finally eluted in low TE buffer.The quality and quantity (region of 150-500 bp) of libraries were assessed using high sensitivity D1000 screen tape on a 4200 TapeStation system.Paired-end sequencing was performed at 2 ng/µL concentration on a NovaSeq 6000 system pooling 25 libraries/S4 flow cell (Illumina, San Diego, CA) with an average of 136 million reads per sample.Raw data have been submitted to NCBI and GSE221921 provides access to all data (https:// www.ncbi.nlm.nih.gov/ gds/?term= GSE22 1921).

Data analysis
FASTQ files corresponding to the forward and reverse reads for 189 samples in total, 96 FM and 93 control were obtained from Illumina BaseSpace and used for analysis.The files were processed using the Trim Galore (Babraham Bioinformatics) and Cutadapt (DOI: https:// doi.org/ 10. 14806/ ej.17.1.200) tools to perform a quality trimming by removing short, low-quality reads and the adapters.RNA-seq reads were mapped to the reference genome (Gencode.v38)and aligned using STAR aligner 22 .Duplicate reads were removed, and uniquely mapped transcripts were selected using Samtools 23 .The TPMs (Transcripts per million) were computed using StringTie 24 corrected by IsoformSwitchAnalyzeR 25 and normalized using TMM (Trimmed Means of M value) 26 .Differential expression analysis was performed using DESeq2 27 , using all annotated genes.Ontological analysis of gene expression was carried out using the Qiagen Ingenuity Pathway Analysis (IPA), and Gene Set Enrichment Analysis (GSEA) 28 .The Interactome analysis was carried out by Pearson Correlation clustering using the 914 most differentially expressed genes for FM1, the 361 most differentially expressed genes for FM2 and the 402 most differentially expressed genes for FM3 and FM4 using Cytoscape 29 , the clusters were determined by AlegroMcode (AllegroViva Corporation, 2011) using default parameters.

Clinical characteristics of FM patients
The median age of patients participating in this study was 48 years of age with overall ages ranging from 28 to 77 years.The median age of the onset of FM was 36 years of age.The median age of control patients was 39 years, with ages ranging from 20 to 69 years.The clinical characteristics of the entire cohort are presented in Table 1.
We performed hierarchical clustering to test whether using clinical characteristics alone, the FM and control cases can be sub-grouped.A near perfect separation of two groups was observed (Fig. 1A) based on clinical symptoms.Only a minority of cases were misclassified: three control subjects: #286, #332, #335 were assigned to the FM group while two FM patients: #028, #078 were assigned to the control group.These five cases were not included in the downstream analysis.
Hierarchical clustering of patient symptoms in the FM cohort indicated that there were three questions that were redundant: (1) patients with poor sleep and insomnia also had memory impairment, (2) depression was related to anxiety and nervousness, and (3) patients who had body pain were more susceptible to have tender points (Fig. 1B).Principal component analysis (PCA) of the symptoms in FM patients indicate that the nonrandom explanation of variance is represented only on the first component and accounts for 18% of the observed variance (Supplementary Fig. 1).This suggests that globally the symptoms of FM patients were homogeneous, except for patient #078 who appeared to be the only outlier marked by the absence of body pain and tender areas but instead reported the presence of depression with physical fatigue.

Cluster analysis of RNA-seq data
Differentially expressed genes were obtained using linear model with the software 'R' (https:// www.R-proje ct.org/).To visualize the results of unsupervised clustering we plotted the logarithm of the TPMs using the heatmap function of 'R' .A total of 1720 differentially expressed transcripts were used to draw the heatmap using the algorithm of TPMs.Of the 90 control cases, 70 formed a tight cluster and 20 were outliers (Fig. 2A, Supplementary Fig. 2).The 20 outliers from the control group dispersed among FM1-3 patients.However, the entire cohort of 94 FM and 90 controls were used for all downstream analyses.The PCA indicated a homogeneous group of 43 patients which we labeled as FM1 and another group of 30 patients labeled as FM2 with unrelated underlying gene expression.The remaining group of 21 patients could be separated into two clusters of 8 and 9 patients each (FM3 and FM4) and 4 outliers (Figs.2B, 6A).For further analysis, we focused on FM1, FM2 and combined FM3 and FM4 as a group.

Analysis of FM1 subgroup
The FM1 subgroup consisted of 43 patients with a similar gene expression pattern.PCA analysis of 480 most DEGs in the 43 FM1 patients and 90 control patients showed a clear separation of patients between FM1 and controls (Fig. 3A).The first component of the PCA encompassed 82% of the total variation indicating that the disease state (Control vs FM1) is the major cause of gene expression difference between these two groups.To understand biological pathways that are specific to the group of 43 FM1 patients, we first performed interactome analysis to identify functional interactions and to pinpoint the DEGs which are most susceptible to being expressed in the same cells.Then we performed IPA on these DEGs to identify the pathways that are associated with these genes.The interactome analysis identified a major cluster composed of 338 DEGs represented in magenta and a smaller cluster of 24 DEGs represented in green (Fig. 3B).The DEGs represented in magenta indicate the presence of a cell (or group of cells) with coordinated gene expression across patients of the FM1 subgroup while www.nature.com/scientificreports/ the DEGs represented in green appear to represent a biological process.We used the genes included in each of these clusters for IPA analysis and found that the significant pathways represented by the major magenta DEGs belonged to extra-cellular matrix genes involved in connective tissue disorders (pulmonary fibrosis, wound healing, cytoskeletal organization, etc.) (Fig. 3C, Table 2).Also, these DEGs pinpoint the presence of upregulated GP6 pathway and the downregulation of Rho GDP-Dissociation Inhibitors (RHODGI) signaling (Fig. 3D).The minor cluster is composed of 24 DEGs that correspond to cell cycle associated genes (Fig. 3E).

Gene set enrichment analysis (GSEA) of FM1 subgroup
Independent of IPA results that searched pathways inside a private database, we also analyzed the whole gene set using GSEA (Gene Set Enrichment Analysis) 28 .This type of analysis ranks the genes from the most differentially expressed to the least differentially expressed and searched among all the known gene sets present in public databases to identify the most enriched gene set.We also analyzed a custom gene set that code for nine proteins composed of the extra-cellular matrix extracted from (Reactome.org):Collagen, Fibrinogen, Elastin, Fibrillin,

Analysis of FM2 subgroup
The FM2 subgroup was composed of 30 patients who showed a similar type of gene expression pattern.Principal component analysis of the 361 DEGs in the 90 control and 30 FM2 patients showed a separation of patients between the controls and FM2 on the first component (Fig. 5A).The first component of the PCA encompassed 34% of the total variation while the second encompassed 6.3% indicating that the disease state (Control vs FM2) is an important cause of gene expression differences between these two groups of patients.Interactome analysis separated the DEGs of the FM2 patients into two distinct clusters (Fig. 5B).These 2 clusters correspond to the upregulated and down-regulated DEGs.IPA analysis of the DEGs identified that the most significant results in this group were the suppression or dysregulation of inflammatory processes (Fig. 5C).The top-ranked dysregulated pathways include phagosome formation, pyroptosis signaling pathway, TREM1 signaling, neuro-inflammation signaling, Th1 pathway, IL1-mediated inhibition of RXR function, crosstalk between dendritic cells and natural   killer cells, toll-like receptor signaling, inflammaosome pathway, Th2 pathway (Fig. 5D, Table 3).These results indicate a lymphocyte to monocyte ratio imbalance in FM2 patients (Fig. 5C).The CLEAR signaling pathway and LXR/RXR activation pathways were upregulated (Fig. 5D).

Analysis of FM3 and FM4 subgroups
The FM3 and FM4 subgroups were composed of 17 cases together.Principal component analysis of the 361 DEGs in the 90 control and 17 FM patients showed a separation of patients between the controls, FM3 and FM4 on the first component (Fig. 6A).The first component of the PCA encompassed 29.5% of the total variation while the second encompassed 4.3% indicating that the disease state (Control vs FM3 and 4) is an important cause of gene expression differences between these groups of patients.Interactome analysis did not establish any significant interactions due to the small sample size.IPA analysis of the DEGs identified the following top-ranked dysregulated pathways including interferon signaling, death receptor signaling, natural killer cell signaling, JAK/STAT signaling and the processing of capped intron-containing pre-mRNA pathway (Fig. 6B, Table 4).

Discussion
FM was previously characterized as a syndrome with widespread pain and localized tenderness 30 .Conceptually, the definition of FM has evolved over time and is perceived as a continuum representing an increased and heightened processing of pain within the nervous system 31 .In 2016, nociplastic pain was proposed as a mechanistic descriptor for FM and chronic pain.Nociplastic pain is defined as pain arising because of an increased sensitivity due to alterations in the peripheral and central nervous system 5 .In FM patients, nocipalstic pain can occur as a comorbidity with an inflammatory, immune, endocrine, genetic and psychosocial factors; all these phenotypes leading to a sensitization phenomenon characterized by a decrease in pain tolerance to afferent nociceptive stimuli 1,32,33 .Over the years, there has been increasing recognition that chronic pain conditions are heterogeneous with a degree of overlap of phenotypes 1,34 .However, to date, there is no clear explanation to account for www.nature.com/scientificreports/this clinical heterogeneity.Our group reported a multiplex cytokine assay that could be used for achieving an objective diagnosis of FM patients 6 .The present study aimed to further define these patients via the identification of genomic markers or signatures to aid diagnosis and possibly lead to the development of mechanism based targeted therapy rather than symptom-based treatment.To this end, we utilized high throughput RNAsequencing for whole transcriptome analysis in 94 patients with FM and 90 healthy control subjects (who had a negative cytokine assay result) using RNA from peripheral blood.The results of our analysis identified multiple subgroups within the cohort of FM patients with distinct non-overlapping gene signatures.Of note, subgroups of FM were identified with enough patients in two subgroups (FM1 and FM2) and combined FM3 and FM4 for detailed downstream analysis.The presence of multiple subgroups within FM patients reflects the inherent clinical heterogeneity associated with FM and chronic pain disorders which explains the diagnostic difficulty often encountered in a clinical setting.The two major subgroups displayed distinct transcriptional profiles indicating two different etiologies that are grouped together under the same general diagnosis of FM.Although we did not identify a specific cause of FM, identification of these subgroups will help develop additional novel diagnostic markers and therapeutics for these patients.The differences observed among the patients suggest that different treatment approaches will be required for patients with FM.
Our study identified subgroups of FM defined by transcriptional signatures.The first group, FM1, included individuals with a signature enriched for gene expression of extracellular matrix (ECM) associated with connective tissue disorders and down regulation of Rho GDP Dissociation Inhibitor (RhoGDI) signaling pathway.The second group, FM2, included individuals that showed a profound reduction in the expression of inflammatory mediators and increased expression of genes involved in the Coordinated Lysosomal Expression And Regulation (CLEAR) signaling pathway.The other two, FM3 and FM4 subgroups, while distinct from the FM1 and FM2, had two few subjects to clearly define the pathways involved.A combined analysis of FM groups 3 and 4 identified overexpression of interferon alpha/beta and JAK/STAT pathways and downregulation of the processing of capped intron containing pre mRNA pathway.In the FM1 subgroup, among the biological processes regulated by the DEGs, include the processing of the ECM protein collagen, wound healing, fibrosis, and genes associated with cell cycle and DNA damage checkpoint regulation (Table 2).Genes associated with the RhoGDI signaling pathway were under-expressed.RhoGDI signaling pathway is the regulator of the Rho family of GTPases that are implicated in the formation of stress www.nature.com/scientificreports/fibers and in pain perception through somatosensory neurons 35 .In this group of patients, deregulation of ECM and tissue homeostasis is likely mediated by fibrocytes.Fibrocytes are bone marrow-derived mesenchymal progenitor cells that directly contribute to tissue remodeling and fibrosis of tissues throughout the body by producing ECM proteins (collagen 1 and collagen III) and by secreting matrix metalloproteinases following injury, wound healing and during fibro-proliferative disorders in response to local tissue injury 36,37 .Fibrocytes traffic to sites of injury during the earliest phase of the innate immune response and exhibit both the inflammatory features of macrophages and the tissue remodeling properties of fibroblasts.They are also an important cellular source of inflammatory cytokines, chemokines, and growth factors that contribute to important autocrine and paracrine signals within the tissue microenvironment 38 .Fibrocytes are distinguished by the simultaneous expression of CD34 or CD45 and collagen 36,37 .Inhibition of Rho kinase increases resting tissue tension which regulates actomyosin contractility, the formation of stress fibers (actin-myosin filaments) and the maturation of focal adhesions 35,39 .Our results suggest that Rho-dependent remodeling of cell matrix is affected in the FM1 subgroup.At the local cellular level, matrix tension has been shown to influence a wide variety of cellular events including neurite growth and angiogenesis 40,41 .Thus, cell-mediated regulation of connective tissue tension may be important to protect blood vessels, sensory and autonomic nerves from prolonged tissue loads induced by various body positions such as sitting, standing, and sleeping positions.In vivo connective tissue tension may not only impact connective tissue homeostasis but also the vascular, nervous, and immune cell populations that reside within the connective tissue network as well as in adjacent organ-specific cell populations.The presence of cell cycle associated genes in the signature indicates persistent stimuli triggered by stress or chronic inflammation that leads to defects in DNA repair mechanisms prompting the activation of fibrocytes 42,43 .
In the FM2 subgroup, there was a significant immune dysregulation as reflected by the under expression of genes involved in phagosome formation, pyroptosis signaling, TREM1 signaling, neuro-inflammation, Th1 and Th2 pathways, crosstalk between dendritic cells and natural killer cells, toll-like receptor signaling and the inflammasome pathway, among others (Table 3).One of the top-ranked pathways that showed overexpression of genes includes the CLEAR signaling pathway.CLEAR pathway is a cellular program that regulates lysosomal biogenesis and participates in macromolecule clearance 44 .CLEAR network is activated by lysosomal storage.The transcription factor EB (TFEB) is a master regulator of lysosomal function 44,45 .TFEB promotes the expression of genes involved in lysosomal biogenesis, such as the mannose 6-phospate receptors, which transport newly synthesized lysosomal enzymes from Golgi to lysosomes.The activity of TFEB is regulated by multiple kinases, in particular the mechanistic target of rapamycin complex 1 (mTORC1) [46][47][48] .When phosphorylated, TFEB is retained in the cytoplasm and inhibited.Several stress signals including nutrient deprivation, proteotoxicity, and lysosomal damage, which have been reported to promote TFEB dephosphorylation, nuclear translocation and activation, leading to an increase in the number and activity of lysosomes 49 .mTORC1 is activated by nutrients and growth factors, and conversely is inhibited by starvation 50 .The activation or inactivation of mTORC1 in response to nutrient availability occurs on the lysosome and is regulated by several lysosomal membrane-associated proteins 51,52 .Thus, the lysosome not only functions as a scaffolding organelle but also participates in the nutrient sensing process.The regulation of mTORC1 signaling by the lysosome also occurs through a transcriptional mechanism mediated by TFEB, which is activated in response to lysosomal stress.TFEB direct target genes were identified by combining ChIP-seq, TFEB overexpression, promoter analysis and coexpression meta-analysis 53 .These genes encode for proteins that can be grouped into several distinct categories, www.nature.com/scientificreports/including lysosomal hydrolases and accessory proteins, lysosomal membrane proteins, subunits of the proton pump, proteins participating in autophagy and non-lysosomal proteins involved in lysosomal biogenesis 53 .Our data show differential expression of genes encoding lysosomal hydrolases and accessory proteins ASAH1, GAA, GNS, IFI30, PSAP; and genes involved in lysosomal acidification ATP6V0B, ATP6V0C, ATP6V0D1, ATP6V1B2 (Table 3) suggesting dysregulation of lysosomal homeostasis in FM2 patients.TFEB also promotes the formation of autophagosomes and their fusion with lysosomes through the upregulation of several key autophagy and lysosomal genes, a process that is initiated by nutrient starvation and executed by the inhibition of extracellular signal regulated kinase 2 (ERK2)-mediated phosphorylation of TFEB at Ser142 48 .Our results show differential expression of the autophagy gene GABA type A receptor-associated protein (GABARAP) (Table 3).GABARAP is a ubiquitin-like modifier that plays a role in intracellular transport of GABA(A) receptors and its interaction with the cytoskeleton.It is involved in autophagy while the microtubule-associated protein 1A/1B-light chain 3 (LC3) is involved in elongation of the phagophore membrane.The GABARAP subfamily is essential for a later stage in autophagosome maturation 54 .Through its interaction with the reticulophagy receptor TEX264, GABARAP participates in the remodeling of subdomains of the endoplasmic reticulum into autophagosomes upon nutrient stress, which then fuse with lysosomes for endoplasmic reticulum turnover 55 .Other TFEB direct targets are genes belonging to distinct families of pattern recognition molecules including membrane-anchored Toll-like receptors (TLRs), which are involved in the innate immune detection of danger signals and microbial motifs 56 and the insulin signaling pathway 53 .Taken together, these results indicate defects in vesicle transport and lysosomal homeostasis in FM2 patients.
In the FM3 and FM4 patients, we identified pathways related to acute inflammatory associated with the Th1 responsive processes with overexpression of interferon pathway, JAK/STAT pathway, IL2, pyroptosis, cell death receptor and necroptosis pathways.Strong down regulation of processing of pre-mRNA pathway indicates global dysregulation of the transcription machinery (Table 4).
There were many limitations in this study inherent to conducting research involving live human subjects and a disease where biology is poorly understood.This study was biased because only individuals who tested positive for the cytokine assay (FM/a) were included.Although this test was developed and validated by our group, it is not widely used or validated by an outside group.Since these criteria were used for patient selection, the study was inherently biased.However, all analyses were performed using the clinical diagnostic criteria for FM recommended by the American College of Rheumatology.We also included only patients with FM rather than patients with chronic pain.In that regard the study is biased towards a subset of patients with FM and the results may not apply to all patients with chronic pain disorder.Future studies are necessary to gain insight into the biological problems in people with chronic pain.
In conclusion, the whole transcriptome analysis of FM patients identified novel gene expression signatures.To our knowledge, this is the first study to report genetic heterogeneity within FM patients.The two major groups of FM patients reported here have defects in the tissue homeostasis associated with ECM and the lysosomal biogenesis pathway.We provide possible mechanisms of FM pathogenesis that need to be further validated to gain precise understanding of the biology of FM and develop novel treatment approaches.

Figure 1 .
Figure 1.Hierarchical clustering of FM patients and controls.(A) Hierarchical clustering of patient symptoms indicates a near perfect separation into two groups.A minority of patients were misclassified: three control subjects #286, #332, #335 were assigned to the FM group while two FM patients #028, #078 were assigned to the control group.Average linkage using Euclidean metrics with k = 2 classes, control patients represented in blue, FM patients represented in red, misclassified patients are identified by an asterisk.(B) Clustering of FM patient symptoms.The symptoms were grouped by hierarchical clustering using average linking of Pearson correlation metric.The vertical bars represent the three groups of correlated symptoms that are present.The star '*' or '0' represents the significance of the cluster.'0': No significant correlation, '*' : adjusted p-value ≤ 0.01.

Figure 2 .
Figure 2. Summary of DEGs identified in FM subgroups.(A) Representation of the proportion of controls and FM patients based on their gene expression.70 control patients with similar gene expression clustered together.The 20 controls that did not have similar gene expression profile are indicated as outliers.Among the FM patients, subgroups with different gene expression profiles were detected.The major group called 'FM1' was formed by 43 patients, a second group called 'FM2' was formed by 30 patients.The rest of the FM patients could be separated into two clusters of 8 and 9 cases respectively and 4 outliers.(B) PCA of FM and controls.PCA of the entire cohort of 94 FM and 90 controls was performed using 1169 DEGs that showed the most significant differential expression.The first component axis shows 46.1% while the second component shows 8.8% of the information.No other components than the first and second components were found useful.Ellipses show 80% confidence interval of each group, the supersized dot corresponds to the centroid of the group.Control patients are represented by a dark blue dot, control patients classified as outliers are represented by a cyan dot, fibromyalgia patients FM1, FM2, FM3, FM4 are represented by a red, green, magenta, orange dot respectively.The 4 FM patients that did not group together are represented by brown dots.

Figure 3 .
Figure 3. Analysis of DEGs in FM1 subgroup.(A) PCA of the 90 controls and the 43 patients of FM1 subgroup using 480 DEGs that showed the most significant differential expression.The near perfect separation shows that most of the variation is represented by the first component (82%).Blue dots: control patients, Red dots: FM patients, the ellipses correspond to the threshold at 80% confidence.Cyan dots show the outlier controls.(B)The interactome analysis of 43 FM1 patients showed a group of 338 DEGs (magenta color) and a small cluster of 24 DEGs (green color).Each dot represents a DEG, the gray dots represent DEGs that did not reach significance.(C) Pathway analysis of the DEGs in the major cluster (magenta) and the minor cluster (green) using a threshold of -log (p value) ≥ 2. (D) Summary of biological functions related to the pathways identified in FM1 patients.(E) Small clusters of 24 genes (green) including coding and lncRNA from the interactome analysis, are associated with cell cycle regulation.Blue: Downregulation, Orange: Upregulation.

Figure 4 .
Figure 4. Results of GSEA analysis of the whole transcriptome for FM1 group.(A) Enrichment of olfactory receptor activity (Normalized Enrichment Score = 2.1, p-value < 0.05).(B) Enrichment of the genes expressing nine proteins associated with the extra cellular matrix (Normalized Enrichment Score = 2.0, p-value < 1E-4).

Figure 5 .
Figure 5. PCA and pathway analysis of FM2 subgroup.(A) PCA of the 90 control and the 30 FM2 patients using 361 DEGs.PCA could separate most of the control and FM2 patients.Blue dots: control patients, green dots: FM2 patients, the ellipses correspond to the threshold at 80% confidence.Cyan dots indicate the outlier controls.(B) Interactome analysis of 361 DEGs that separated the FM2 using the 90 control patients as reference into two distinct clusters of up and down regulated DEGs.Blue: down regulated genes, Red: up regulated genes.(C) IPA analysis shows the most significant pathways found in FM2 until -log (pvalue) ≥ 3. The gene expression indicates numerous pathways associated with inflammatory response are downregulated and the CLEAR signaling pathway is upregulated.(D) Summary of biological functions related to the pathways identified in FM2 patients.Blue: downregulation, Orange: upregulation.

Figure 6 .
Figure 6.PCA and pathway analysis of FM3 and FM4 subgroups.(A) First component axis is showing 29.5% and was used for separation of control patients from FM3 and FM4.Ellipses show 80% confidence interval of each group, and the supersized dots correspond to the centroid of the group.Control patients are represented by a dark blue dot, control outliers are represented by a cyan dot, FM3 are represented by a magenta dot while FM 4 are represented by an orange dot.The 4 outlier FM cases were represented by brown dots.(B) IPA analysis shows the most significant pathways found in FM3&4 until -log (pvalue) ≥ 2. The gene expression indicates numerous pathways associated with acute inflammatory processes are upregulated and the processing of pre-mRNA pathway is downregulated.

Table 1 .
Clinical characteristics of FM and control groups.

Table 2 .
Top ranked canonical pathways in FM1.

Table 3 .
Top ranked canonical pathways in FM2.

Table 4 .
Top ranked canonical pathways in FM3 and FM4.